#==============================================================================#
# 6-georef.R
# Matt Curtis mjdcurtis@gmail.com
# updated Feb 14, 2025
# updated May 29 2025
#------------------------------------------------------------------------------#
# old geocoding method
# later we switched to photon
#==============================================================================#

rm(list = ls())
gc()

library(stringr)
library(arrow)
library(tidytable)
library(progress)
library(RecordLinkage)

#library(tictoc)
#------------------------------------------------------------------------------#
rm(list = ls())
gc()
#tic()

# setwd once
setwd("~/Replication package")
genidir<-"~/Replication package/"


#------------------------------------------------------------------------------#
# apply the same string cleaning to each location

cleanloc <- function(loc) {
  loc <- loc %>%
    iconv(to = 'ASCII//TRANSLIT') %>% # remove accents
    str_to_lower() %>%
    # remove everything but letters, - and _
    # then replace - with _
    str_replace_all("[^A-za-z_\\- ]", " ") %>% 
    str_replace_all("\\[", " ") %>% 
    str_replace_all("\\]", " ") %>% 
    # str_replace_all("\\-", "_") %>%
    str_squish() %>%
    replace(is.na(.), "")
  return(loc)
}


#------------------------------------------------------------------------------#
# loading geonames text dump

geonames<-fread("./1_raw_data/1_2_geni/auxiliary_files/FR.txt", 
                     sep = "\t") %>% tidytable

names(geonames)<-c(
  'geonameid',
  'name',
  'asciiname',
  'alternatename',
  'latitude',
  'longitude',
  'feature class',
  'feature_code',
  'country code',
  'cc2',
  'admin1_code',
  'admin2_code',
  'admin3_code',
  'admin4_code',
  'population',
  'elevation',
  'dem',
  'timezone',
  'modification_date'
)

#------------------------------------------------------------------------------#
# geonames administrative regions
# all regions above a city

admin1<-geonames %>% 
  filter(feature_code=="ADM1") %>% 
  select(name,admin_code=admin1_code) %>% 
  mutate(admin_level = 1)

altadmin1<-geonames %>% 
  filter(feature_code=="ADM1") %>% 
  select(alternatename,admin_code=admin1_code) %>% 
  separate(alternatename,into="name",sep=",")%>% 
  mutate(admin_level = 1)

admin2<-geonames %>% 
  filter(feature_code=="ADM2") %>% 
  select(name,admin_code=admin2_code)%>% 
  mutate(admin_level = 2)

altadmin2<-geonames %>% 
  filter(feature_code=="ADM2") %>% 
  select(alternatename,admin_code=admin1_code) %>% 
  separate(alternatename,into="name",sep=",")%>% 
  mutate(admin_level = 2)

admin3<-geonames %>% 
  filter(feature_code=="ADM3") %>% 
  select(name,admin_code=admin3_code)%>% 
  mutate(admin_level = 3)

altadmin3<-geonames %>% 
  filter(feature_code=="ADM3") %>% 
  select(alternatename,admin_code=admin1_code) %>% 
  separate(alternatename,into="name",sep=",")%>% 
  mutate(admin_level = 3)

admin4<-geonames %>% 
  filter(feature_code=="ADM4") %>% 
  select(name,admin_code=admin4_code)%>% 
  mutate(admin_level = 4)


altadmin4<-geonames %>% 
  filter(feature_code=="ADM4") %>% 
  select(alternatename,admin_code=admin1_code) %>% 
  separate(alternatename,into="name",sep=",")%>% 
  mutate(admin_level = 4)


geonames_admin<-bind_rows(admin1,admin2,admin3,admin4) %>% 
  bind_rows(altadmin1,altadmin2,altadmin3,altadmin4)
rm(admin1,admin2,admin3,admin4)
rm(altadmin1,altadmin2,altadmin3,altadmin4)

# remove some helper words
for(remove in c("Arrondissement","Departement","de",
                "du","d'","de l'","des")){
  geonames_admin <-geonames_admin %>% 
    mutate(name = str_replace_all(name,remove,"")) %>% 
    mutate(name = str_squish(name))
}

geonames_admin<-mutate(geonames_admin,name = cleanloc(name)) 
gc()

#------------------------------------------------------------------------------#
# geonames cities

geonames_cities<-geonames %>% 
  filter(`feature class`=="P") %>% 
  mutate(city = cleanloc(name)) %>% 
  select(latitude,longitude,city)%>% 
  filter(city!="")%>% 
  mutate(city = str_replace_all(city,"_"," "))

geonames_alt<-geonames %>% 
  # mutate(geoid = row_number()) %>% 
  # select(alternativename,geoid) %>% 
  separate(alternatename,into="alternative",sep=",")%>% 
  filter(`feature class`=="P") %>% 
  mutate(city = cleanloc(alternative)) %>% 
  select(latitude,longitude,city) %>% 
  filter(city!="")

geonames_cities<-bind_rows(geonames_cities,geonames_alt) %>% distinct() %>% 
  mutate(city = str_replace_all(city,"_"," "))

#------------------------------------------------------------------------------#
# geni data

output <-
  read_parquet(paste0(genidir, '2_scripts/2_0_tempfiles/fr-raw.parquet')) %>% 
  tidytable() 

output$rowid<-1:nrow(output)


# I had not previously merged in marriage place

marriages<-read_parquet(paste0(genidir,'2_scripts/2_0_tempfiles/fr-pids-m.parquet')) %>% 
  tidytable() %>% 
  select(married_france,profile_id) %>% 
  distinct() 

before<-nrow(marriages)
output<-left_join(output,marriages)
stopifnot(nrow(marriages)==before)

#------------------------------------------------------------------------------#
# location variables

# births
fr_births <- output %>%
  filter(born_france==TRUE) %>% 
  select(birth_city,
         birth_county,
         birth_state,
         birth_country,
         birth_place_name
  ) %>% 
  distinct() %>% 
  mutate(
    birth_city = cleanloc(birth_city),
    birth_county = cleanloc(birth_county),
    birth_state = cleanloc(birth_state),
    birth_country = cleanloc(birth_country),
    birth_place_name = cleanloc(birth_place_name),
  ) %>% 
  distinct() 

names(fr_births)<-c("city","county","state","country","place")

# deaths
fr_deaths <- output %>%
  filter(died_france==TRUE) %>% 
  select(birth_city,
         birth_county,
         birth_state,
         birth_country,
         birth_place_name
  ) %>% 
  distinct() %>% 
  mutate(
    birth_city = cleanloc(birth_city),
    birth_county = cleanloc(birth_county),
    birth_state = cleanloc(birth_state),
    birth_country = cleanloc(birth_country),
    birth_place_name = cleanloc(birth_place_name),
  ) %>% 
  distinct() 

names(fr_deaths)<-c("city","county","state","country","place")

# deaths
fr_marriages <- output %>%
  filter(married_france==TRUE) %>% 
  select(birth_city,
         birth_county,
         birth_state,
         birth_country,
         birth_place_name
  ) %>% 
  distinct() %>% 
  mutate(
    birth_city = cleanloc(birth_city),
    birth_county = cleanloc(birth_county),
    birth_state = cleanloc(birth_state),
    birth_country = cleanloc(birth_country),
    birth_place_name = cleanloc(birth_place_name),
  ) %>% 
  distinct() 

names(fr_marriages)<-c("city","county","state","country","place")

fr_locs <- fr_births %>%
  bind_rows(fr_deaths) %>% 
  bind_rows(fr_marriages) %>% 
  #select(loc) %>% 
  distinct() 

fr_locs$locid<-1:nrow(fr_locs)

fr_locs$locid<-1:nrow(fr_locs)

nloc <- nrow(fr_locs)

# use this later
write_parquet(fr_locs,"./2_scripts/2_0_tempfiles/fr-locations-original.parquet")

#------------------------------------------------------------------------------#
# donation step
# if there is no city, take from county etc.
#------------------------------------------------------------------------------#

# donate cities
# define a potential city as anything with 2+ letters 
print("donate city")

cities <-fr_locs$city %>%
  unique() %>% 
  tidytable() %>%
  filter(. != "") %>% 
  mutate(len=str_length(.)) %>% 
  filter(len>1)

 
fr_nocity<-filter(fr_locs,city=="") %>% 
  # I use these later
  mutate(oldflag=999) %>% 
  mutate(match=NA) %>% 
  # mush together the other variables
  mutate(mush = paste(place,county,state,sep=" | "))

# search for one of the cities from the possible list in mush
# if more than one, take the first
pb <- progress_bar$new(total = nrow(cities)+1)
for(c in cities$.){
  # followed by end of line or space or pipe
  regex<-paste0(c,"(?=($|\\|| ))")
  # look for location of regex in mush
  fr_nocity$flag<-str_locate(fr_nocity$mush,regex)[,1]
  # replace with 999 if NA
  
  # if the newly found potential city comes earlier in mush, replace match and
  # oldflag 
  fr_nocity$flag<-replace(fr_nocity$flag,is.na(fr_nocity$flag),999)
  fr_nocity<-fr_nocity %>% 
    mutate(match = case_when(
      flag<oldflag~c,
      TRUE ~ match
    ))%>% 
    mutate(oldflag = case_when(
      flag<oldflag~flag,
      TRUE ~ oldflag
    ))
  pb$tick()
}
rm(pb)
gc()


fr_nocity<-fr_nocity %>% 
  # clean up by pulling out the components of mush and removing the extracted match
  # I am not sure if this is 100% working as intended
  mutate(place = str_extract(mush,"^[^|]*(?=\\|)")  %>% 
                    str_remove(match) %>% str_squish) %>% 
  mutate(county = str_extract(mush,"(?<=\\|)[^|]*(?=\\|)") %>% 
           str_remove(match) %>% str_squish)%>% 
  mutate(state = str_extract(mush,"(?<=\\|)[^|]*$") %>% 
           str_remove(match) %>% str_squish) %>% 
  mutate(city=match) %>% 
  select(-flag,-oldflag,-match,-mush)%>% 
  mutate(city_donated=1)


# reassemble and check no duplicates
before<-nrow(fr_locs)
fr_locs<-filter(fr_locs,city!="") %>% 
  mutate(city_donated=0) %>% 
  bind_rows(fr_nocity)
stopifnot(nrow(fr_locs)==before)

#------------------------------------------------------------------------------#
# donate county
print("donate county")

fr_locs<-mutate(fr_locs,county=replace(county,is.na(county),""))

# define a potential city as anything with 2+ letters 
counties <-fr_locs$county %>% unique() %>% 
  tidytable() %>%
  filter(. != "") %>% 
  mutate(len=str_length(.)) %>% 
  filter(len>1)

fr_nocounty<-filter(fr_locs,county=="") %>% 
  # I use these later
  mutate(oldflag=999) %>% 
  mutate(match=NA) %>% 
  # mush together the other variables
  mutate(mush = paste(place,state,sep=" | "))
  
# search for one of the counties from the possible list in mush
# if more than one, take the first
pb <- progress_bar$new(total = nrow(counties)+1)
for(c in counties$.){
  # followed by end of line or space or pipe
  regex<-paste0(c,"(?=($|\\|| ))")
  # look for location of regex in mush
  fr_nocounty$flag<-str_locate(fr_nocounty$mush,regex)[,1]
  # replace with 999 if NA
  fr_nocounty$flag<-replace(fr_nocounty$flag,is.na(fr_nocounty$flag),999)
  # if the newly found potential coounties comes earlier in mush, replace match and
  # oldflag 
  fr_nocounty<-fr_nocounty %>% 
    mutate(match = case_when(
      flag<oldflag~c,
      TRUE ~ match
    ))%>% 
    mutate(oldflag = case_when(
      flag<oldflag~flag,
      TRUE ~ oldflag
    ))
  pb$tick()
  
}
pb$tick()
rm(pb)
gc()

# clean up by pulling out the components of mush and removing the extracted match
# I am not sure if this is 100% working as intended
fr_nocounty<-fr_nocounty %>% 
  mutate(place = str_extract(mush,"^[^|]*(?=\\|)")  %>% 
           str_remove(match) %>% str_squish) %>% 
  mutate(state = str_extract(mush,"(?<=\\|)[^|]*$") %>% 
           str_remove(match) %>% str_squish) %>% 
  mutate(county=match) %>% 
  select(-flag,-oldflag,-match,-mush)%>% 
  mutate(county_donated=1)


# reassemble and check no duplicates
before<-nrow(fr_locs)
fr_locs<-filter(fr_locs,county!="") %>% 
  mutate(county_donated=0) %>% 
  bind_rows(fr_nocounty)
stopifnot(nrow(fr_locs)==before)

#------------------------------------------------------------------------------#
# donate state
print("donate state")


fr_locs<-mutate(fr_locs,state=replace(state,is.na(state),""))

states <-fr_locs$state %>% unique() %>% 
  tidytable() %>%
  filter(. != "") %>% 
  mutate(len=str_length(.)) %>% 
  filter(len>1)

  
fr_nostate<-filter(fr_locs,state=="") %>% 
  # I use these later
  mutate(oldflag=999) %>% 
  mutate(match=NA) %>% 
  # mush together the other variables
  mutate(mush = paste(place,state,sep=" | "))


# search for one of the counties from the possible list in mush
# if more than one, take the first

pb <- progress_bar$new(total = nrow(states)+1)
for(c in states$.){
  regex<-paste0(c,"(?=($|\\|| ))")
  fr_nostate$flag<-str_locate(fr_nostate$mush,regex)[,1]
  fr_nostate$flag<-replace(fr_nostate$flag,is.na(fr_nostate$flag),999)
  # if the newly found potential coounties comes earlier in mush, replace match and
  # oldflag 
  fr_nostate<-fr_nostate %>% 
    mutate(match = case_when(
      flag<oldflag~c,
      TRUE ~ match
    ))%>% 
    mutate(oldflag = case_when(
      flag<oldflag~flag,
      TRUE ~ oldflag
    ))
  pb$tick()
  
}
rm(pb)
gc()
print("DONE")

# clean up by pulling out the components of mush and removing the extracted match
# I am not sure if this is 100% working as intended
fr_nostate<-fr_nostate %>% 
  mutate(place =mush %>% 
           str_remove(match) %>% str_squish) %>% 
  mutate(state=match) %>% 
  select(-flag,-oldflag,-match,-mush)%>% 
  mutate(state_donated=1)

# reassemble and check no duplicates
before<-nrow(fr_locs)
fr_locs<-filter(fr_locs,state!="") %>% 
  mutate(state_donated=0) %>% 
  bind_rows(fr_nostate)
stopifnot(nrow(fr_locs)==before)

#------------------------------------------------------------------------------#
# link to admin codes
#------------------------------------------------------------------------------#
# here I take a detour to link
# all the possible city, county, states to admin codes
# the goal is to deal with a specific centroid bias
# the one from being born in "antwerp, france" when that is 
# an arrondissement not a city

# all the cities
cities <-fr_locs$city %>% unique() %>% 
  tidytable() %>%
  filter(. != "") 
names(cities)<-c("name")
cities$type<-"city"

# all the counties
counties <-fr_locs$county %>% unique() %>% 
  tidytable() %>%
  filter(. != "") 
names(counties)<-c("name")
counties$type<-"county"

# all the states
states <-fr_locs$state %>% unique() %>% 
  tidytable() %>%
  filter(. != "") 
names(states)<-c("name")
states$type<-"state"

# all the places combined
locs<-bind_rows(cities,counties,states) %>% distinct()
locs$id<-seq(1:nrow(locs))


#------------------------------------------------------------------------------#
# 1.a.
# link on direct match

direct<-left_join(locs,geonames_admin) %>% 
  mutate(matches = n(),.by=c(id))

badzero<-filter(direct,is.na(admin_level)) %>% 
  select(name,id,type) %>% distinct()

goodtemp<-filter(direct,!is.na(admin_level)) %>% 
  select(name,type,id,admin_code,admin_level)%>% distinct()

good<-goodtemp %>% select(id,type) %>% distinct()


print(paste0("1+: ",nrow(good)))
print(paste0("0: ",nrow(badzero)))
print(paste0("Gained ",nrow(good)+nrow(badzero)-nrow(locs)," obs."))

#------------------------------------------------------------------------------#
# 1.b
# fuzzy match
print(nrow(badzero))
fuzzy_out<-data.table()

# do a fuzzy string match using jarro winkler with a threshold of .9
pb <- progress_bar$new(total = nrow(badzero)+1)
for(i in 1:nrow(badzero)){
  name<-badzero$name[i]
  id<-badzero$id[i]
  dist<-jarowinkler(name,geonames_admin$name) %>% tidytable
  names(dist)<-"jwscore"
  
  dist<-mutate(dist,score= case_when(
    jwscore>=.9 ~ 1,
    jwscore<.9 ~ 0
  ))
  
  geonames_admin$jwscore<-dist$jwscore
  out<-geonames_admin[geonames_admin$jwscore>=.9,]
  out$old_name<-name
  out$id<-id
  out$matches<-nrow(out)
  
  fuzzy_out<-bind_rows(fuzzy_out,out)
  pb$tick()
  
}
pb$tick()
rm(pb)
gc()


good_fuzzy<-filter(fuzzy_out,!is.na(admin_level))
goodtemp<-bind_rows(goodtemp,good_fuzzy)

badzero<-left_join(badzero,fuzzy_out) %>% filter(is.na(matches))


good<-goodtemp %>% select(id,type) %>% distinct()

print(paste0("1+: ",nrow(good)))
print(paste0("0: ",nrow(badzero)))
# I don't think this is a problem but it does gain some here
print(paste0("Gained ",nrow(good)+nrow(badzero)-nrow(locs)," obs."))

#------------------------------------------------------------------------------#
# find admin codes
# 

goodtemp<-select(goodtemp,-type) %>% distinct() 

goodtemp<-goodtemp%>% 
  filter(admin_level!=4)

before<-nrow(fr_locs)

# look in state for an admin div. name
state_admin <- fr_locs %>%
  left_join(select(
    goodtemp,
    state = name,
    state_level = admin_level,
    state_code = admin_code
  )) %>%
  arrange(locid, state_level) %>%
  mutate(state_n = row_number(), .by = locid) %>%
  filter(state_n == 1) %>% 
  select(state_level, state_code, locid)

fr_locs<-left_join(fr_locs,state_admin)

# look in county for an admin div. name
# but drop it if there is a similar div. in the county
county_admin <- fr_locs %>%
  left_join(select(
    goodtemp,
    county = name,
    county_level = admin_level,
    county_code = admin_code
  )) %>%
  filter((county_level != state_level) | is.na(state_level)) %>%
  arrange(locid, county_level) %>%
  mutate(county_n = row_number(), .by = locid) %>%
  filter(county_n == 1) %>% select(county_level, county_code, locid)

fr_locs<-left_join(fr_locs,county_admin)

# look in city for an admin div. name
# but drop it if there is a similar div. in the state or county
city_admin <- fr_locs %>%
  left_join(select(
    goodtemp,
    city = name,
    city_level = admin_level,
    city_code = admin_code
  )) %>%
  filter((city_level != county_level) | is.na(county_level)) %>%
  filter((city_level != state_level) | is.na(state_level)) %>%
  arrange(locid, county_level) %>%
  mutate(county_n = row_number(), .by = locid) %>%
  filter(county_n == 1) %>% select(city_level, city_code, locid)

fr_locs <- left_join(fr_locs, city_admin)


stopifnot(nrow(fr_locs)==before)

# for each level, pick the most local where it shows up
fr_locs<-fr_locs %>% 
  mutate(admin_code3=case_when(
    city_level ==3 ~ city_code,
    county_level ==3 ~ county_code,
    state_level ==3 ~ state_code,
    TRUE ~ ""
  ))%>% 
  mutate(admin_code2=case_when(
    city_level ==2 ~ city_code,
    county_level ==2 ~ county_code,
    state_level ==2 ~ state_code,
    TRUE ~ ""
  ))%>% 
  mutate(admin_code1=case_when(
    city_level ==1 ~ city_code,
    county_level ==1 ~ county_code,
    state_level ==1 ~ state_code,
    TRUE ~ ""
  ))

#------------------------------------------------------------------------------#
# 1. link!
#------------------------------------------------------------------------------#


cities <-fr_locs$city %>% unique() %>% 
  tidytable() %>%
  filter(. != "") 

nocities<-fr_locs %>% filter(is.na(city))
print(paste0(nrow(nocities),"/",nrow(fr_locs)," have missing cities and will not be matched"))

names(cities)<-c("city")
cities<-cities %>% 
  mutate(city = str_replace_all(city,"_"," "))

cities$cityid<-seq(1:nrow(cities))


#------------------------------------------------------------------------------#
# 1.a.
# link on direct match

direct<-left_join(cities,geonames_cities) %>% 
  mutate(matches = n(),.by=c(cityid))

# two or more matches is bad
bad2plus_temp<-filter(direct,matches>1) 

bad2plus<-bad2plus_temp%>% 
  select(cityid,city) %>% distinct()

# zero matches is bad
badzero<-filter(direct,is.na(longitude)) %>% 
  select(cityid,city) %>% distinct()

# one match is good!
good<-filter(direct,!(matches>1|is.na(longitude))) %>% 
  select(cityid,city,latitude,longitude)%>% distinct()

print(paste0("1: ",nrow(good)))
print(paste0("2+: ",nrow(bad2plus)))
print(paste0("0: ",nrow(badzero)))
print(paste0("Gained ",nrow(bad2plus)+nrow(good)+nrow(badzero)-nrow(cities)," obs."))

#------------------------------------------------------------------------------#
# 1.b
# fuzzy match using a jw score of 0.9 

print(nrow(badzero))
fuzzy_out<-data.table()

pb <- progress_bar$new(total = nrow(badzero)+1)
for(i in 1:nrow(badzero)){
  city<-badzero$city[i]
  cityid<-badzero$cityid[i]
  dist<-jarowinkler(city,geonames_cities$city) %>% tidytable
  names(dist)<-"jwscore"
  
  geonames_cities$jwscore<-dist$jwscore
  out<-geonames_cities[geonames_cities$jwscore>=.9,]
  out$city<-city
  out$cityid<-cityid
  out$matches<-nrow(out)
  
  
  fuzzy_out<-bind_rows(fuzzy_out,out)
  pb$tick()
  
}
pb$tick()
rm(pb)
gc()

# reassemble the three categories, good, bad with 2+, and bad with 0
# note that the case with 2+ plus fuzzy matches but 1 direct match,
# we will not have done the fuzzy match to begin with 
good_fuzzy<-filter(fuzzy_out,!is.na(latitude)&matches==1)
good<-bind_rows(good,good_fuzzy)

badzero<-left_join(badzero,good_fuzzy) %>% filter(is.na(matches))
bad2plus_fuzzy<-filter(fuzzy_out,!is.na(latitude)&matches>1)
bad2plus_temp<-left_join(bad2plus_temp,bad2plus_fuzzy)

bad2plus<-bad2plus_temp%>% 
  select(cityid,city) %>% distinct()

print("")
print(paste0("1: ",nrow(good)))
print(paste0("2+: ",nrow(bad2plus)))
print(paste0("0: ",nrow(badzero)))
print(paste0("Gained ",nrow(bad2plus)+nrow(good)+nrow(badzero)-nrow(cities)," obs."))


#------------------------------------------------------------------------------#
# 2. for two plus, look at the admin codes!

bad2plus_temp <- left_join(bad2plus_temp, geonames) %>%
  select(
    city,
    cityid,
    latitude,
    longitude,
    geoadmin1 = admin1_code,
    geoadmin2 = admin2_code,
    geoadmin3 = admin3_code
  )

fr_locs_good<-fr_locs %>% left_join(good) %>% filter(!is.na(latitude))

fr_locs_2plus<-fr_locs %>% left_join(bad2plus_temp) %>% 
  filter(!is.na(latitude)) %>% 
  filter(geoadmin1==admin_code1|admin_code1=="") %>% 
  filter(geoadmin2==admin_code2|admin_code2=="") %>% 
  filter(geoadmin3==admin_code3|admin_code3=="") %>% 
  # keep the ones we manage to elimiate this way
  mutate(dups = n(),.by=locid) %>% 
  filter(dups==1)

fr_locs_good<-fr_locs_good %>% 
  bind_rows(fr_locs_2plus) %>% 
  mutate(flag = n(),.by=locid)

stopifnot(max(fr_locs_good$flag)==1)
# how many did we geolocate?
print(nrow(fr_locs_good)/nrow(fr_locs))


# export with quality control flags
final<-read_parquet("./2_scripts/2_0_tempfiles/fr-locations-original.parquet") %>% 
  tidytable() %>% 
  left_join(select(fr_locs_good,locid,city_donated,latitude,longitude,city_level)) %>% 
  mutate(city_admin_flag = as.numeric(!is.na(city_level))) %>% 
  mutate(city_donated_flag = as.numeric(city_donated==1))

# use this later
write_parquet(final,"./2_scripts/2_0_tempfiles/fr-geocoded-v3.parquet")
#toc() %>% print()